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Abstract 

The direct numerical simulation (DNS) of the Taylor-Couette flow in the fully tur- 
bulent regime is described. The numerical method extends the work by Quadrio & 
Luchini (Eur. J. Mech. B / Fluids, 21, 413-427, 2002), and is based on a parallel 
computer code which uses mixed spatial discretization (spectral schemes in the ho- 
mogeneous directions, and fourth-order, compact explicit finite-difference schemes in 
the radial direction) . A DNS is carried out to simulate for the first time the turbulent 
Taylor-Couette flow in the turbulent regime. Statistical quantities are computed to 
complement the existing experimental information, with a view to compare it to 
planar, pressure-driven turbulent flow at the same value of the Reynolds number. 
The main source for differences in flow statistics between plane and curved-wall 
flows is attributed to the presence of large-scale rotating structures generated by 
curvature effects. 



1 Introduction 



The flow in tlie gap between coaxial rotating cylinders, known as the Taylor- 
Couette flow (TC flow hereinafter), is among the most investigated problems in 
fluid mechanics (see [1] for an overview), owing to its engineering applications 
[2], as well as to its relevance as prototypical flow in the study of transition 
to turbulence and of fully-developed turbulent flows over streamwise-curved 
surfaces. 

Many experimental and numerical papers have dealt with the instabilities 
developed by the TC flow when the value of the Reynolds number Re is 
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increased such that the laminar regime is taken over by a sequence of bifur- 
cations, eventually leading to chaotic behavior. Less is known about the fully 
turbulent regime, since a smaller number of laboratory experiments and es- 
sentially no direct numerical simulation (DNS) studies are available. The lack 
of DNS investigations can be ascribed - at least in part - to the difficulties 
of implementing an efficient method for the numerical solution of the incom- 
pressible Navier-Stokes equations in cylindrical geometries, and at the same 
time to the computational demand of such simulations in turbulent regime. 
Hence a satisfactory quantitative description of the turbulent TC flow, lever- 
aging the full space-time information potentially available from a DNS, is as 
yet missing. 

In the first part of this paper a numerical method for the DNS of the Navier- 
Stokes in cylindrical coordinates is presented, which is designed for computing 
the Taylor-Couette flow in the turbulent regime. The method rehes on the 
strategy developed by Quadrio & Luchini [3] for the DNS of a turbulent flow 
in an annular pipe: the main enhancements introduced here are an improved 
accuracy of the spatial discretization, and the ability of the code to exploit 
parallel computing on commodity hardware. 

The newly developed numerical method is then used to carry out what we be- 
lieve is the first DNS of the turbulent Taylor-Couette fiow. We shall focus on a 
so-called small-gap geometry, identical to that considered by Andereck, Liu & 
Swinney [4] , and set the Reynolds number to a relatively high value. Analysis 
of the results will focus on the assessment of low-order flow statistics, and on 
some theoretical issues concerning turbulence in presence of wall curvature. As 
described at length in the review paper [5] by Patel & Sotiropoulos, curvature 
considerably impacts the modeling of turbulent flows. Its effects thus carry 
a broad interest, and also affect experimental practices. For example, when 
measuring the skin friction in a turbulent flow over a curved wall, one often 
resorts to methods (like the Clauser plot) which imply the validity of the law 
of the wall, and which often require the numerical values of its parameters (i.e. 
the slope of the logarithmic part of the profile, given by the inverse of the von 
Karman constant, and its intercept) to be known in advance. Indeed, neither 
these values possess to date an agreed-upon value, nor is their dependence on 
the degree of curvature known. A further example occurs in the derivation of 
the friction law for the TC fiow, i.e. a formula relating the friction coefficient 
to the value of the Reynolds number: such derivation, actively discussed in 
literature [6,7,8,9], is often based on assumptions (logarithmic form vs. power 
law) about the shape of the mean velocity profile. A DNS, with its ability to 
evaluate the skin friction directly, is a powerful tool to complement experimen- 
tal data in this research areas. Of course the main limitation of DNS lies in 
the limited values of the Reynolds number typically affordable in the numeri- 
cal simulations. This is the reason why in this paper the numerical method is 
developed with great emphasis on computational efficiency. 
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A DNS at a relatively high value of Re will also offer the opportunity to ob- 
serve to what extent a planar, pressure-driven and fully-developed turbulent 
flow (like a channel flow) presents analogies to the TC flow, in the spirit of 
what has been already observed by Lee & Kim [10] for the planar TC flow. 
The present simulation reaches an unprecedently high value of Re, so that 
the comparison is made at the value Rer ~ 180 (based on the friction ve- 
locity and half the gap width), which is considered the minimum bound for 
the near-wall viscous turbulent cycle to be fully active. A significant physi- 
cal process which acts in the TC flow and is absent in plane geometries is 
the additional mass and momentum transfer exerted by large-scale toroidal 
rolls. The laminar solution is known to be unstable above a well-defined (and 
curvature-dependent) critical value Rcc of the Reynolds number, above which 
such structures, known as Taylor vortices, are quickly generated [11]. As Re 
increases further, the Taylor vortices undergo a series of transformations, after 
which they eventually reappear in the turbulent regime [12]. In his illuminated 
paper [13], Townsend surmised that turbulence in the TC flow is basically of 
two different kinds, with one contribution from the wall shear and the other 
from the large-scale structures. DNS is the perfect tool to analyze the full flow 
fleld and to describe these large-scale structures, that are absent both in the 
plane channel and in the plane Couette flows. 

The outline of the paper is as follows. In §2 the geometry of the problem is 
introduced, and the main elements of its numerical simulation are recalled 
in §2.1 (additional information is deferred to an Appendix). A validation of 
the computer code and a critical comparison of its results to available data 
will be given in §3, together with a quantitative assessment of the computa- 
tional performance. The numerical simulation of the turbulent TC flow will 
be described in §4, with regard to discretization parameters and computa- 
tional procedures. In §5 results will be presented in terms of averaged and 
instantaneous flow properties. Lastly, §6 will contain a conclusive summary. 



2 Problem definition and numerical method 

We consider the classic Taylor-Couette apparatus, made by a moving inner 
cylinder with radius TZi rotating at an angular velocity Qi, and a concentric, 
outer fixed cylinder with radius TZo- In figure 1 the geometry employed in 
the present work is sketched: the radial coordinate is r, the axial and az- 
imuthal coordinates being x and 9. The corresponding velocity components 
are v (wall- normal) , u (spanwise) and w (streamwise) . The amount of cur- 
vature is expressed by the geometric parameter ( = TZi/TZo- The Reynolds 
number Re is defined with the inner cylinder rotating speed Wi = VL{IZi, the 
gap width IZo — TZi = 2h between the two cylinders and the kinematic vis- 
cosity y of the fluid. For a given value of C, when gravity and other external 
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Fig. 1. Geometry of the Taylor-Couette flow system. The inner cyUnder rotates at 
speed Wi and angular velocity fij, the outer cylinder is at rest. 

forces are neglected, Re becomes the sole relevant fluid dynamic parameter. 
Unless otherwise stated, the results presented in the following are made non- 
dimensional by using Wi as reference velocity and h as reference length. When 
useful, the wall distance y is also used, conveniently defined as y = r — TZi oi 
y — TZo — r. Quantities made non-dimensional in wall units are indicated with 
a -|- superscript and are computed with friction velocity Ur (see §5 for details 
concerning its definition) and fiuid viscosity u. 



2.1 Governing equations and discretization 

The Taylor-Couette system is of course most easily described in cylindrical 
coordinates. On the other hand, a widely used and extremely efficient formu- 
lation of the Navier-Stokes equations exists in cartesian coordinates, and is 
well suited to the DNS of turbulent fiows with two homogeneous directions. 
This formulation was introduced by Kim, Moin & Moser [14] , and is based on 
rewriting the equations of motion in terms of two scalar equations for the wall- 
normal components of the velocity and vorticity vector. The computational 
efficiency is maximized when flow variables are expanded in Fourier series in 
the two homogeneous directions (the wall-normal discretization is irrelevant). 



4 



A few studies have described the extension of this formulation to a cyhn- 
drical geometry, to enjoy the same advantages: pressure removed from the 
equations, and optimal computational load. We build for the present work 
upon the contribution by Quadrio & Luchini [3], who wrote evolutive equa- 
tions for the radial components of velocity and vorticity in cylindrical coordi- 
nates, and used finite differences to compute radial derivatives with schemes 
of second-order formal accuracy. The time integration was carried out in [3] 
with an explicit third-order Runge-Kutta method for the nonlinear terms and 
curvature-related viscous terms, and a second-order Crank-Nicholson scheme 
for the remaining viscous terms. 

The present work employes a numerical method which is improved in at least 
two respects compared to Ref. [3]: the second-order accuracy of the finite- 
difference schemes was not considered satisfactory, and - most importantly - 
distributed-memory parallel computing capabilities are required. As shown by 
Luchini and Quadrio in their recent paper [15], these two issues can be dealt 
with in an unified approach. The formal accuracy of the radial derivatives - 
and, most importantly, their spectral resolution - is improved by resorting to 
high-accuracy compact difference schemes. The schemes are moreover made 
explicit, by following the procedure illustrated by Thomas [16] which leverages 
the absence of the third derivative in the equations of motions, to reduce 
their computational cost. At the same time, finite differences, which enjoy the 
property of being local operators in physical space, are key to obtain a good 
parallel performance on commodity networking hardware. 

Reproducing the approach of Ref. [15] in the cylindrical case, however, is not 
trivial. As shown in the Appendix, the governing equations require further 
manipulation in order to arrive at a form where explicit compact schemes 
can be applied. A computational stencil made by 5 arbitrarily spaced (and 
smoothly stretched) grid points is used to obtain a formal accuracy of order at 
least 4, and the coefficients are computed, according to an approach common 
in the theory of Pade approximants, by matching the discrete derivative of a 
function known in analytic form (e.g. a polynomial) to its analytical derivative. 
We refer the reader to Ref. [15] for a general illustration of the method in the 
cartesian case, and to the Appendix for a short discussion of the peculiarities 
of the cylindrical case. Full description of the method is given by Ref. [17]. 



3 Validation 

The present computer code is carefully validated against both numerical and 
experimental data available for the TC fiow through preliminary calculations 
at non-turbulent values of Re. 
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Fig. 2. Left: computed laminar velocity profile (symbols) and laminar analytical 
solution (1) (continuous line), for ^ = 0.5 and A^^ = 33 points in the radial direction. 
Right: difference between computed and analytical solutions, as a function 

of Nr. The straight line shows the expected decrease proportional to N~'^. 

The correctness and accuracy of finite-difference radial operators is preliminar- 
ily checked by computing the exact laminar solution. In the laminar regime, 
the flow is described by an exact solution of the Navier-Stokes equations [18]: 



Due to the term ~ r~^, this solution cannot be represented exactly by the 
polynomial interpolation implied by a FD method. A discretization error is 
thus always present, which is bound to decrease as the step size raised to the 
fourth power. The laminar velocity profile computed for ( — 0.5 is shown in 
Fig. 2 together with the analytical solution. The difference between computed 
and analytical curves is indeed observed to decrease as requested by the formal 
fourth-order accuracy of the numerical method. 

The entire code is then tested both in large- and small-gap geometries for 
the first insurgence of Taylor vortices, to verify that the evolution of small- 
amplitude disturbances is well represented. The curvature-dependent critical 
value of the Reynolds number, above which instability amplifies fiow per- 
turbations to develop Taylor vortices, is computed with the present code and 
compared to available data. For a given C, the code is run at various val- 
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ues of Re, starting from an initial condition made by the laminar sohition 
with superimposed small-amplitude random disturbances. The discretization 
parameters are reported in table 1 ( "straight- vortex" case) . The temporal evo- 
lution of the kinetic energy is monitored, to verify whether the initial energy 
decreases or gets amphfied depending on Re. For ( — 0.5 we find Rec to be 
confined within the range 68.2 < Rcc < 68.4. At the same curvature, previ- 
ous numerical simulations by Fasel & Booz [19], based on an axisymmetrical 
finite-difference method, determined Rec = 68.2, whereas experimental mea- 
surements [20] indicated Rec — 68.4. For the small-gap case, at C = 0.95 we 
compute 184 < Rec < 186, which agrees with the value of 185 determined by 
Moser et al. [21]. 

In the same flow regime, the torque needed to maintain the rotation of the 
inner cylinder is computed as an indicator of global spatial accuracy, and 
compared with previous works. The torque is defined as: 



where L is the axial length of the cylinders. 

Table 2 show how values of T computed by the present code at C = 0.5 
compare very well (within small fractions of a percent) to literature data. 

A further increase in Re brings us into the non-linear regime, where the 
straight Taylor vortices undergo a wavy azimuthal deformation and are named 
"wavy vortices" after the paper [22] by Coles. The prediction of quantitative 
properties of the wavy Taylor vortices, like their azimuthal wavenumber and 
their phase speed, allows a further confirmation of the present code. A simu- 
lation with the discretization parameters reported in table 1 ( "wavy- vortex" 
case) has been carried out at ^ = 0.5 and Re = 250. In this simulation only, to 
ease the comparison with the paper [23] by Snider, the outer cylinder possesses 
a small amount of counter-rotation, such that Re based on the outer cylinder's 
rotation speed is Re — —55. The simulation, described in detail in [17], shows 
a pair of wavy vortices, with azimuthal wavenumber m = 2, and a charac- 
teristic rotation period T = IVJh/Wi. This translates into a non-dimensional 
rotation speed s — ^-K/mTVti — 0.054, which compares very well to the results 



3. 1 Computational performance 

Here the computational performance of the code is documented for the prob- 
lem size and on the computing system used for the turbulent simulation de- 
scribed later in §4. 





of Ref. [23]. 
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We employ a spatial resolution of A^^ = 512 and A^^ = 170 Fourier modes in 
the azimuthal and axial directions respectively, while the number of points in 
the radial direction is set to Nr = 129. Parallel simulations are carried out on a 
dedicated computing system, called a Personal Supercomputer, which is based 
on commodity personal computers. The special architecture of this system, 
that can be used with either the cartesian or the cylindrical code, has been 
described in Ref. [15]. The system is made by 10 SMP nodes, each equipped 
with 2 Intel Xeon 2.66 GHz CPU, and at least 1GB of 266 MHz SDRAM. The 
nodes carry two Gigabit Ethernet adapters each, and are connected in a ring 
topology. 

On one Xeon CPU the code requires 680 MB of memory, and takes 182 seconds 
to compute a full time step for a 3-substeps Runge-Kutta temporal scheme. 
The amount of memory includes storage space for the radial coefficients dis- 
cussed in the Appendix. Should RAM size become an issue, these coefficients 
can can be computed on-the-fly, thus saving 13% of RAM at the expense of a 
comparable increase in computing time. The CPU overhead of the cylindrical 
code compared to its cartesian counterpart is about 40% overall. 

Table 3 documents how the computation is speed up by employing an increas- 
ing number A^ of computing nodes. The ratio between the wall-clock time on 
one CPU and the wall clock time with A^ distributed-memory CPUs, i.e. the 
parallel speedup S, is not far from the linear one: the PLS parallel strategy 
in the present context thus exploits fully our low-cost computing system. The 
maximum measured speedup with A^ = 10 is about 9. The SMP capabili- 
ties of the nodes are then exploited on top of the distributed-memory parallel 
strategy. This allows a further 1.6 speedup factor from the use of 2 CPUs. 
The SMP speedup is additive to the distributed-memory speedup, and shows 
essentially no decrease with increasing the number of computing nodes. The 
entire wall clock time for the simulation described in §4 can thus be decreased 
from about 6 months to about 11 days when the entire system is used. 



4 The turbulent Taylor— Couette flow simulations 

The geometry considered in the present analysis is the small-gap geometry 
described by Andereck, Liu & Swinney [4] with ( = 0.882, which gives TZi = 
15h and TZo = I7h. The value of the Reynolds number is set at Re = 10500, to 
obtain a value of the friction Reynolds number similar to Re-r — 180, which is 
typically considered the lowest value at which a turbulent channel flow presents 
a well-developed inner layer [14,24]. To our knowledge, this is the highest value 
of Re ever reached in a DNS of Couette flow, and of course represents a big 
computational challenge: for comparison, Bech et al. [25] computed a fully 
turbulent plane Couette flow, which is the hmit for C — > 1 of the present TC 
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flow, at Rcr ~ 80 only. According to Hamilton, Kim & Waleffe [26], the plane 
Couette flow is turbulent above Rct — 30. 

The size of the computational domain must be chosen by keeping in mind that 
periodic boundary conditions are used in the two homogeneous directions: a 
periodic box imposes an artificial large- wavelength cutoff to the structures that 
can be represented in the numerical simulation. It can be recalled however that 
periodicity in the streamwise direction is a physically sound condition for the 
TC flow, unhke planar channel and Couette flows. The azimuthal extension 
of the box is chosen so that Lg — IOtt/i when measured at the center line of 
the gap. The turbulent TC flow presents a dominant periodic pattern in the 
axial direction too: toroidal structures present an axial wavelength of 5/i, and 
experimental results [12] indicate that this wavelength is rather insensitive 
to the value of Re. The simulations are thus carried out for — 5h. It 
can be mentioned moreover that the implied indeflnite mirroring of periodic 
boxes makes the present simulation totally free from the end-effects that are 
unavoidable in laboratory experiments. 

To obtain the resolution of all the signiflcant spatial scales, Ng = 512 and 
— 170 Fourier modes are used in the azimuthal and axial directions re- 
spectively, and for the radial direction A^^ = 129. The spatial resolution in 
wall units matches or exceeds the typical resolution employed in channel-flow 
DNS [24]. The global number of degrees of freedom amounts to ^ 2.2 • 10^. 
The parameters related to the spatial discretization are reported in table 4. 
The resolution of the relevant temporal scales dictates the time step size. We 
use At = 0.012S /Wi, which is smaller than the stability limit of the employed 
Rungc-Kutta scheme. In viscous units, this corresponds to At"*" ^ 0.08. A null 
mean pressure gradient is imposed in both the homogeneous directions. 

The adequacy of the present size for the computational domain can be flrst 
evaluated in terms of inner (viscous) units, i.e. by making lengths non-dimensional 
with the friction velocity Ur (see later §5 for its deflnition) and the fluid viscos- 
ity. With the friction velocity computed at the inner wall = 946, whereas 
at the outer wall = 840. The streamwise length is between = 5576 at 
the inner wall and — 5600 at the outer wall. A comparison with the do- 
main sizes usually employed for DNS in the turbulent plane channel flow show 
that the domain size is certainly adequate for representing the inner turbulent 
layer. The streamwise length is 2.5 times the streamwise length typically 
employed [14] in the DNS of a plane turbulent channel flow at a comparable 
value of friction Reynolds number; the spanwise width is large enough to rep- 
resent well the underlying turbulent flow, being larger than that used in Ref. 
[24]. 

It remains to be determined whether representing a single pair of periodic 
structures is actually enough to obtain domain-independent low-order statis- 
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tics. Given the emphasis of the present work towards high possible Re, we 
decided to use the available computational resources to reach the highest Re 
while representing a single pair of TV, and to check a posteriori that low- 
order statistics are not significantly affected by the axial domain size. The 
check consisted in running an additional simulation for a domain size with 
doubled L^, and an accordingly doubled N^. At the same time, the averaging 
time is halved to maintain the computational load approximately constant. 
In short, the result of this check, that will be further discussed in the follow- 
ing, is that representing a single pair of structures is adequate for measuring 
low-order statistics. 



^ . 1 Computational procedures 



The initial condition for the simulations is based on the laminar solution (1); 
divergence-free velocity disturbances with amplitude O(10~^) and random 
phase are added to the whole set of Fourier modes. Different initial condi- 
tions have been considered in (less resolved) preliminary simulations, by either 
changing the amplitude of the disturbances, or by applying disturbances to a 
partial set of modes, or by starting the run from a null mean velocity profile 
instead of the profile given by Eq. (1). In terms of the long-time mean turbu- 
lent friction and others higher-order statistics, no difference has been noticed 
among these cases after the turbulent regime sets up, even though the time 
required to reach the statistically steady-state, as well as the behaviour of the 
fiow during the initial transient, have been observed to depend on the details 
of the initial conditions. 



The simulation is started from the above-described initial condition and let 
run for 10005/1^1. A flow fleld is stored on disk every 10 time units for later 
postprocessing. The size of a single flow fleld is 185 MBytes. 



The time required for the simulation to settle to a statistically steady-state 
requires an initial time interval: we have estimated it based on the time his- 
tory of various quantities, including the derivative of the space-mean velocity 
longitudinal velocity profile friction computed at the walls, which is related to 
Rer- Figure 3 shows how the length of this transient is about 150 time units. 
In computing statistical quantities, fiow fields corresponding to the initial 300 
time units are rejected. 
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Fig. 3. Behaviour of Rct at the inner wall (continuous line) and outer wall (dashed 
line) as a function of time after the startup of the simulation. 

5 Results 

The first global flow feature that we discuss is the mean friction. The shear- 
stress is computed for both the inner and the outer walls. Its mean value 
(r-u,), averaged over time and over the homogeneous directions, allows us to 
compute a friction velocity Uj- = \J (tw) / p local to each wall. From Ur a local 
Reynolds number Rcr = UrS/v can be obtained. The present results yield 
Rtr — 189.3 for the inner wall and Rtr — 167.7 for the outer wall. As a further 
confirmation of the adequacy of spatial resolution and statistical averaging, 
we observe that the ratio between the values of Re^ at the two walls equals 
C, within 0.18%. The check with larger yielded values of Rcr which differ 
from the above ones by less than 1%, i.e. less by the error implied by the finite 
averaging time. 

The radial profile {w) of the mean azimuthal component of the velocity vector 
is shown in figure 4 and compared to the laminar solution (1). While the 
laminar profile at this curvature level is little different from a straight line, 
the turbulent profile presents, in analogy to the plane turbulent Couette flow, 
two distinct regions: near the walls a shear-driven boundary layer, and in 
the central part of the gap a region where the velocity decreases very slowly 
with r. In particular, the large central region presents an almost constant or 
very slowly-increasing angular momentum r{w)] its value of nearly 0.57liWi 
agrees with the measurements by Taylor [11] and Smith & Townsend [27]. 
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Fig. 4. Radial distribution of the mean azimuthal velocity {w) (closed symbols, 
continuous line) and of the angular momentum r{w) (open symbols, dashed line). 
Lines correspond to the analytical laminar solution. Note that angular momentum 
is made non-dimensional by WilZi. The virtually indistinguishable line drawn along 
the closed circles is the {w) profile for the simulation with = lO/i. 

This property of the flow, observed here for the first time in a numerical 
simulation, qualitatively supports the hypothesis of a core region of constant 
angular momentum flow put forward by Townsend [28]. The figure confirms 
moreover that the mean profile is essentially unaffected by the number of 
represented TV pairs: the profile computed for the simulation with — lOh 
is hardly distinguishable from the one with — 5h. 

In figure 5 the azimuthal mean velocity profile is plotted in semi-logarithmic 
scale and in local wall units against the wall distance y'^. Similarly to pressure- 
driven flows, a viscous sublayer can be observed to develop over both walls for 
^ 5, where the profile follows the linear law (w)^ = widely accepted 
for pressure-driven flows. The profiles over both walls are fairly similar up to 
y+ fa 40 (buffer layer), but in the central part of the gap the outer profile 
becomes slightly higher than the inner one: this indicates that the friction 
velocity (local to each wall) alone is not the correct scaling velocity in this 
region of the flow. 

The very existence of a logarithmic layer, where the mean velocity profile is 
described by: 

{w)+ = -lny+ + S, 
as well as its characterization through the numerical values of the von Karman 
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Fig. 5. Mean azimuthal velocity profile {w)~^ , in local wall units, over the inner 
wall (closed symbols) and the outer wall (open symbols). The continuous line is the 
linear profile (w)'^ = 

constant k, and the intercept B, are matters scholars have not yet agreed upon 
[5]. Figure 6 shows the azimuthal mean velocity profile, plotted in law-of-the- 
wall form for the inner wall only. The computed profile is compared to a 
few laws proposed in the literature. In particular for a plane Couette flow at 
Rcr = 82.2 the values 1/k, = 2.55 and B = 4.7 have been suggested [25], 
whereas for the same flow at Rcr = 52 the values 1/k, = 2.5 and B = 4.6 have 
been used [29]. The values for the plane channel flow at Rct = 180, after Kim, 
Moin & Moser [14], are 1/k = 2.5 and B = 5.5. 

Figure 6 reveals a rather limited extent of the logarithmic region, the relatively 
high value of Re notwithstanding. (A rough estimate is that the present proflle 
has constant slope for 20 < |/+ < 40. The outer edge of the logarithmic layer 
can also be estimated from figure 5 as the position where the curves relative to 
the two walls start to diverge. On the contrary, the profile of the plane channel 
flow has constant slope in a range of nearly 100 wall units.) The values of 1/k 
and B from the plane channel case do not provide a good fit, whereas the ones 
from the plane Couette flow appear more suited to the present velocity proflle. 
The best fitting line seems to be that by Ref . [25] . The main information that 
can be drawn is thus to reinforce the statement by Smith & Townsend [27]: 
"no significant region of logarithmic variation of velocity can exist" for "any 
flow of Reynolds number less than 20000" in a Taylor-Couette flow with the 
outer cylinder at rest. According to a semi-empirical law [30], which relates 
the outer-scale Reynolds number to Rct, this corresponds to Rcr 340. As 
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Fig. 6. Mean azimuthal velocity profile {w)~^ (symbols) over the inner wall, compared 
with logarithmic laws proposed in literature. Dashed line: plane Couette flow at 
Rer = 52 by [25]. Dash-dotted line: plane Couette flow at Rcr = 82.2 by [29]. 
Continuous line: plane channel flow at Rer = 180 by [24]. 

long as such information is not available, the collected evidence suggests that, 
at least for these relatively low values oi Re, k should take values equal or 
slightly reduced compared to the plane channel flow. We recall however that 
curvature for the considered geometry at C = 0.882 is rather weak. 

The observed form of the mean velocity profile is relevant to the ongoing 
discussion about the derivation of a friction law for the turbulent Taylor- 
Couette flow. In recent years, Lathrop et al [6] derived an approximate friction 
law based on the assumption of a logarithmic velocity profile, whereas Panton 
[7] started from the assumption that the core region has constant angular 
momentum. The present results, which show to a very good approximation 
the presence of the constant angular momentum region, appear to lend more 
support to the latter assumption. However, tough wider logarithmic layers 
could develop at higher i?e, our data do not rule out the possibility that 
different theories (see for example Ref. [9] and [31]) may yield a better friction 
law. 

The variance of the velocity fluctuations is plotted in figure 7 across the whole 
gap. One observes the highest levels of fiuctuations for the streamwise com- 
ponent, which is also the one with more marked asymmetry between the two 
walls and shows higher turbulence activity over the outer wall. The radial com- 
ponent peaks at the centerline with a symmetric profile, whereas the spanwise 
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Fig. 7. Variance of velocity fluctuations across the gap. Triangles: streamwise compo- 
nent. Diamonds: wall-normal component. Circles: spanwise component. For better 
clarity only one every two points is shown. 

component is nearly symmetric and remarkably intense. It emerges clearly 
how the wall-normal component completely differs from the other ones in the 
central region of the gap. 

A different look can be obtained by plotting the root-mean-square values for 
both walls in local wall units against in figure 8, thus implicitly accounting 
for the different values of the friction Reynolds number. Convex and concave 
walls do still present discernible differences: in particular the turbulence ac- 
tivity as deduced by the level of streamwise velocity fiuctuations is still larger 
over the outer wall. A similar but less pronounced behaviour is observed for 
the spanwise component. Inner scaling seems to work for the radial compo- 
nent, at least very near the wall: when the centerline is approached, the two 
profiles over the inner and outer wall appear to diverge. 

We interpret the failure of friction velocity in scaling the r.m.s. profiles and 
the particular behaviour of the v component as effects of the turbulent Taylor 
vortices. Near the wall they induce perturbations in the u and w components 
without affecting v, owing to the presence of the solid wall. In a sense, these 
perturbations behave similarly to the so-called inactive motions hypothesized 
by Townsend [28] to populate the turbulent boundary layer, and known to 
be responsible for the so-called anomalous scaling [32,33]. Figure 8 reports 
also the r.m.s. value of velocity fluctuations for the plane channel flow, taken 
from [24]. The streamwise component is quite similar, while the wall-normal 
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Fig. 8. Root-mean-square value of velocity fluctuations above both walls. Closed 
symbols: inner wall. Open symbols: outer wall. Lines without symbols: plane channel 
flow at Rcr = 180 by [14]. Triangles and continuous line: streamwise component. 
Diamond and dashed line: wall-normal component. Circles and dotted line: spanwise 
component. For clarity only one every two points is shown. 

component is similar in the near-wall region only. The spanwise component 
in the channel flow presents a much lower level of activity, and is thus the 
component where the contribution by the turbulent Taylor vortices is the 
most relevant. 

Figure 9 reports a snapshot of the whole flow fleld. A pair of turbulent Taylor 
vortices can be clearly observed. In laboratory experiments, they are usually 
visualized by resorting to passive tracers, whereas DNS makes all the flow 
variables easily available. Fig.9 shows the entire computational domain: two 
isosurfaces where the radial velocity component attains a value of ±0.12VFj 
suggest the presence of two large-scale roll-like structures. It is evident how the 
vortices are strongly modulated by the noisy turbulent background, so that 
their boundary is somewhat blurred. They are however clearly identifiable, 
and extend down to the wall. 



6 Discussion and conclusions 

A direct numerical simulations of a Taylor-Couette flow in the fully turbulent 
regime has been presented. The value of the Reynolds number - up to Rcr ~ 
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Fig. 9. Snapshot of an instantaneous flow field. The bottom waU moves from left to 
right. Isosurfaces for the radial velocity component are shown, at levels of ib0.12Wj: 
negative levels are in light gray, and positive in dark gray (blue). 

190 based on friction velocity, viscosity of the fluid and half gap width - is 
high enough to compare turbulent statistics of the present flow and of planar 
pressure- driven flows. 

Such a demanding simulation has required developing a numerical method de- 
signed ad hoc. It discretizes the governing equations with high accuracy and 
exploits the computing power of a parallel computing system, which is assem- 
bled from commodity hardware. The code, based on Fourier expansion in the 
homogeneous directions and on fourth-order, compact flnite differences in the 
radial direction, has been vahdated by computing several physical quantities 
available from previous DNS or experiments in the laminar and transitional 
regimes, and has then made possible the flrst DNS of the Taylor-Couette flow 
in the fully turbulent regime. 

A geometry with relatively small curvature has been considered. In comparing 
statistics to those from other wall turbulent flows, analogies and differences 
have been observed, the main source for the latter being the presence of the 
curvature- induced turbulent Taylor vortices. The mean velocity proflle does 
not exhibit an equilibrium logarithmic region at the present value of Re, con- 
trasting the pressure-driven flows. A core region with almost constant angular 
momentum has been revealed by our simulations, and the numerical value of 
this constant compares very well to experimental observations taken at differ- 
ent Re. The analysis of the root-mean-square values of velocity fluctuations, 
in comparison with the ones from a plane channel flow, reveals the lack of 
inner scaling, thus pointing to the existence of different physical processes 
contributing to turbulent fluctuations. 
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Two distinct sources of fluctuations can indeed be identified in the present flow: 
(i) the large-scale vortices, generated by an instability mechanism ultimately 
related to curvature, and (ii) the wall turbulence cycle, linked to the near- 
wall shear, which in the present context produces small-scale fluctuations and 
provides a noisy background for the vortices. A visualization of a snapshots of 
the flow field has shown how the large-scale structures are embedded in the fiow 
and retain their basic shape, though blurred by the surrounding small-scale 
turbulence. The interaction of the small-scale turbulence with the large-scale 
structures, and the possibility of discerning and eventually separating their 
contributions to the flow statistics is an interesting long-term objective on 
which we are currently working. 
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Appendix 

This Appendix provides some details concerning the entire procedure of ex- 
tending the numerical method from what documented in Ref. [3] to the one 
employed in the present simulations. The complete derivation can be found in 



Since two directions are homogeneous, the first step consists in transforming 
the governing equations in Fourier space. Fourier-transformed quantities will 
be denoted with a hat; a and m denote the axial and azimuthal wavenum- 

bers, and k'^ = a'^ + /r'^. Note that k depends on the radial coordinate. 
In [3], after manipulation of the Navier-Stokes equations in primitive vari- 
ables the following equation for the Fourier-transformed radial vorticity fj was 
determined: 



Ref. [17]. 



dfj 



1 




dt 



Re 
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In Eq. (3) the so-called Chandrasekhar notation is employed for radial deriva- 
tives: 



D{f) 



dr ' 



dr r ' 



and HU, HW and (in the following Eq.(4)) HV denote the non-hnear terms 
as written in the original primitive- variables equations. 

After further manipulation, and leveraging the continuity equation, the tempo- 
ral evolution of v was shown in [3] to be described by the following fourth-order 
differential equation: 



d_ 

at 



Re 



D^DD^{v) - 2—v + 2—D{w) - 2—w 



+ 



1 

Re 



D 



1 



— tm — \ 
ia HU + —HWj 



+ HV. (4) 



The entire differential system requires 6 boundary conditions; the no-slip and 
no-penetration conditions imposed at both walls imply v = 0, D{v) — and 
fi — at r — TZi and r — TZg [34] . 

Unlike the cartesian case, it can be noted that u and w appear in Eqs. (3) 
and (4) through some of the viscous terms too (note that such terms vanish 
in the planar limit r — > oo), owing to the form of the Laplacian operator in 
cylindrical coordinates. Hence, when a partially-implicit procedure is used for 
temporal integration, these curvature-related viscous terms cannot be solved 
for implicitly. The solution adopted in [3] was to solve for them explicitly, 
since this was shown to imply no stability limitations. 

Though a FD discretization is more convenient than a spectral one in terms of 
parallel computing, one has to deal with its inferior spatial accuracy. The ac- 
curacy can be made comparable to that of spectral methods by using compact, 
high-accuracy FD schemes [35] . A further advantage is that the computational 
cost of compact schemes, that are usually implicit and require solving a system 
to evaluate a derivative, can be further reduced for the Navicr-Stokes equa- 
tions: in the present context, they can be made explicit, as observed as early 
as in 1953 by Thomas [16], when third derivatives are absent in the governing 
equations. We thus begin by eliminating third derivatives in Eq. (4) via the 
continuity equation. The simple substitution: 

7^ /-\ - im ^ 
D^[v) — —lau w. 
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does the job, and of course causes a few viscous terms to be moved into the 
part of the evohition equation (4) that is to be integrated exphcitly in time. 
Similarly to the previously mentioned terms, this does not limit temporal 
stabihty, and the choice of the time step size is still dictated by considerations 
of temporal accuracy. 

Now, if a{r) indicates a generic r-dependent coefficient and / is the unknown 
function, terms in the form aD[f) do not lend themselves to a straightforward 
use of compact explicit schemes as in [15]. Such terms must be rewritten after 
repeated integrations by parts, i.e. with substitutions of the type: 

aD{f) = D{af) - D{a)f; aD,{f) = D,{af) - D{a)f. (5) 



Note that integrating by parts introduces additional known coefficients like 
D{a), depending on the radial coordinate only (directly and/or via the wavenum- 
bers): one of the simplest among them is D(l/k'^). 

The procedure of repeated integrations by parts needed to massage Eqns. (3) 
and (4) into a form suitable for the application of the compact finite-difference 
operators is now described in some detail. 

In Eqn. (4), the first term which is integrated by parts following (5) is the 
time derivative, which becomes: 



In the right-hand-side of Eqn. (4), perhaps the most complicated term is: 

■ 1 



D 



^2 [-D,DD,v) 



where first the continuity equation must be invoked to cancel the third deriva- 
tive. Repeated integrations by parts then allow the r-dependent coefficients 
to remain only in the innermost positions. After some algebra, the result is: 



D 



^ {D,DD,v) 



-DD^DD, [^v]+D 



1 



-DD —V 



- 2DD. 



-D 



+ D 



DDD 



V 



1 

SDD 



r 

-D 
D 



1 



:D 



+ 



+ 



im 

lau H w 

r J 



Having used the continuity equation, the last term has appeared which cannot 
be solved for implicitly, and must thus join the viscous curvature terms in the 
explicit part. 
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The last term of Eqn. (4) which needs manipulation is: 
1 



D 



2im < DD 



w 



-D 



W 



The same sequence of integration by parts must be carried out for Eqn. (3) 
for the radial vorticity, arriving at the following substitution: 



tTTl 

2 — D(u) = 2im 
r 



D[-]+2- 



Lastly, the nonlinear terms too contain radial derivatives, and some terms 
therein must be integrated by parts. 

This procedure leads to the final, rather long form of the equations for v and 
r), which lends itself to a discretization in the radial direction with explicit 
compact finite difference schemes of fourth-order accuracy over a five point 
stencil. It is written here in full for completeness, without time discretization 
for notational simplicity: 



d_ 

di 

±\ 
Re I 

+D 



v-DD, 



2DD^{v) - DD^DD^ 
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+ 2im 



DD^(-uv)-Di^D y-^uv 



DD 
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+3D 
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VW 



+D 



1 / ^am 



DD 

uw 4 
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+ 
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\d(^]w 



fc2 



+ 



m 



-w 



+ 



lauv 



D{v'^) - —vw - -v^ + -w^; (6) 



dfi _ 1 
dt ~ Re 



DD^{fi) - k'^fj + 2im 



Di—u)+2—-—v 



— im 



la 



^ ( uv\ 2 im 

D \ — + —uv H — -uw 



+ ia 



lauw + D \ vw) H w H — vw 

r r 



■ (7) 



At this stage, the strategy for efficient parallel computing that has been suc- 
cessfully used in our cartesian Navier-Stokes solver, recently described as the 
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Pipelined Linear System (PLS) method by Luchini and Quadrio in [15], can be 
implemented quite easily. The resulting code can run without modifications 
on the computing machine designed and built for our cartesian code, with 
obvious advantages. 
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Le iV^ Nr Ne Re 
"straight-vortex" case Ah nh 64 64 32 60-190 

"wavy-vortex" case 4/i 27r/i 64 64 32 255 

Table 1 

Size of the computational domain and grid resolution for the preliminary simulations 
in the "straight-vortex" and "wavy- vortex" regimes (see text). 



Re 


here 


C = 0.5 [19] 


60.0 


16.7551 


16.7551 


68.0 


16.7551 


16.7551 


70.0 


17.1542 


17.1537 


75.0 


18.1634 


18.1627 


80.0 


19.0536 


19.0527 



Table 2 

Comparison between the value of the torque T computed by the present code and 
Ref. [19] for C = 0.5. 



s 


1 


1.97 2.8 


3.7 


4.6 


5.5 


6.1 


7 7.8 


8.5 


N 


1 


2 3 


4 


5 


6 


7 


8 9 


10 



Table 3 

Parallel speedup S of the code versus the number N of the computing nodes. 



Position 




Le 


Lt 




Az+ 




(rAO) 


inner wall 


5h 


9.377r/i 


946 


5576 


5.6 


0.9 


10.9 


centreline 


5h 


IOtt^ 


889 


5588 


5.2 


4.6 


10.9 


outer wall 


5h 


10.627r/i 


839 


5600 


4.9 


0.8 


10.9 



Table 4 

Size of the computational domain and grid resolution for the turbulent Taylor- 
Couette simulation. Wall units are computed on the basis of the friction velocity 
relative to each wall, as defined in §5. 
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